Reionization, SLOAN, and WMAP: is the Picture Consistent? 
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ABSTRACT 

I show that advanced simulations of cosmological reionization are able to fit the observed data 
on the mean transmitted flux in the hydrogen Lyman-alpha line at z ~ 6. At the same time, 
posteriori models can be constructed that also produce a large value (20%) for the Thompson 
scattering optical depth, consistent with the WMAP measurements. Thus, it appears that a 
consistent picture emerges in which early reionization (as suggested by WMAP) is complete by 
z ~ 6 in accord with the SLOAN data. 



Subject headings: cosmology: theory 
- galaxies: intergalactic medium 

Introduction 



cosmology: large-scale structure of universe - galaxies: formation 



Are we knocking at the door of reionization? 
At least, it seems so, as the SLOAN quasar sur- 
vey continues to find ever more distant quasars 
(Becker et al. 2001; Fan et a. 2002, 2003; White et 
al. 2003). In fact, highest redshift quasars show es- 
sentially no transmitted flux just blueward of the 
quasar Hn region, which has been interpreted as 
the evidence that that reionization took place at 
a redshift z ~ 6 (Becker et al. 2001; White et al. 
2003; a similar claim was also made by Djorgovski 
et al. 2001 based on an extrapolation from a lower 
redshift observations). 

However, one should be cautious before drawing 
such a conclusion. Indeed, the observed decrease 
in the mean transmitted flux at z ~ 6 might sim- 
ply indicate a decrease in the mean ionizing inten- 
sity rather then a real reionization of the universe. 
After all, a neutral fraction of only 2 x 10~ 4 is suf- 
ficient to absorb effectively all Lyman-alpha radia- 
tion at z = 6 even in the gas that is underdense by 
a factor of 10. Thus, without an understanding of 
the evolution of the universe around the reioniza- 
tion epoch, the Ly-alpha absorption data cannot 
be used to constrain the epoch of reionization (un- 
less damping wings are observed in the absorption 
profiles, Miralda-Escude 1998). 

Fortunately, our theoretical understanding of 



the process of reionization, in part based on nu- 
merical simulations, is solid enough so that the 
evolution of the ionizing intensity at these red- 
shifts can be predicted with a reasonable confi- 
dence level. Combining simulations with the ob- 
servational data indeed allows one to come up with 
a consistent picture for the process of cosmological 
reionization. 

2. Simulations 

Four sets of simulations have been performed 
with the SLH code and are similar to the simu- 
lations reported in Gnedin (2000). The main dif- 
ference with previous simulations is that a newly 
developed and highly accurate Optically Thin Ed- 
dington Variable Tensor (OTVET) approximation 
for modeling radiative transfer (Gnedin & Abel 
2001) is used instead of a crude Local Optical 
Depth approximation. The new simulations there- 
fore should be sufficiently accurate (subject to 
the usual limitations of numerical convergence and 
phcnomenological description of star formation) to 
be used meaningfully in comparing with the ob- 
servational data. In particular, my fiducial model 
used the cosmological parameters as determined 
by the WMAP satellite (Spergel, et al. 2003). 

Parameters of the four sets are given in Table 1 . 
All simulations included 128 3 dark matter parti- 
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cles, an equal number of baryonic cells on a quasi- 
Lagrangian moving mesh, and about 3 million stel- 
lar particles that formed continuously during the 
simulation. The nominal spatial resolution of sim- 
ulations with the box size of 4ft -1 Mpc was fixed 
at lh^ 1 comoving kpc, with the real resolution 
being a factor of two worse. Simulations with the 
box size of 8h^ x Mpc had the twice worse spatial 
resolution. 

In all cases a flat cosmology was assumed, with 
^A,o = 1 — Qm.o, an d normalization of the primor- 
dial fluctuations was determined either from the 
WMAP data (Spergel et al. 2003) for sets A4 and 
A8, or from the COBE data (White & Bunn 1995). 
Notice that a small change in the slope of the pri- 
mordial power-law spectrum n makes a significant 
effect on the amount of the small-scale power due 
to a large leverage arm from CMB scales to the 
tens-of-kpc scales which are important for reion- 
ization. 

Star formation is incorporated in the simu- 
lations using a phenomenological Schmidt law, 
which introduces two free parameters: the star 
formation efficiency esF (as defined by eq. (1) of 
Gnedin 2000) and the ionizing radiation efficiency 
euv (defined as the energy in ionizing photons per 
unit of the rest energy of stellar particles). 

The star formation efficiency esF is chosen so 
as to normalize the global star formation rate in 
the simulation at z = 4 to the observed value from 
Steidel et al. (2001), whereas the ultraviolet radi- 
ation efficiency euv is only weakly constrained by 
the (highly uncertain) mean photoionization rate 
at z ~ 4. The redshift of reionization strongly 
depends on euv and is not in fact predicted in 
a simulation, but can be changed over a reason- 
able range depending on the assumed value of euv- 
Each of the simulations sets included in Table 1 in 
fact included several individual simulations with 
different values of euv- 

However, since the simulations are quite expen- 



Table 1: Simulation Parameters 



Set 


^ m 


h 


n 


Box size (h 1 Mpc) 


A4 


0.27 


0.71 


1.0 


4 


A8 


0.27 


0.71 


1.0 


8 


B4 


0.35 


0.70 


0.95 


4 


C4 


0.35 


0.70 


0.97 


4 
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Fig. 1. — Mean free path to ionizing radiation 
(bottom) and its time derivative (top) as a func- 
tion of redshift for the simulation A4 from Table 
1. 

sive, it is not possible to cover a large range of euv 
in a given set. Typically, only 2 or 3 simulations 
per set have been performed and the results are 
then interpolated between the simulations. This 
procedure is fully described in §4. 

3. A "Redshift of Reionization" : What Is 
R? 

Reionization is a process, not an event. In fact, 
whole process of reionization is quite extended 
(Az ~ 5 — 10), and can be generically separated 
into three stages: (i) the pre-overlap stage, in 
which individual Hn regions around the sources 
of ionization expand and merge in the low density 
IGM, (ii) the overlap stage, in which all indi- 
vidual Hn regions overlap, and last remnants of 
the neutral low density gas quickly disappear, and 
(iii) the post-overlap stage, in which remaining 
high density gas is being ionized from the outside, 
until neutral gas remains only in some of the high- 
est density regions, which would be identified as 
Lyman-limit system in the absorption spectra of 
distant quasars. 

One can think reionization as complete, when 
the mean free path to the ionizing radiation is fully 
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Fig. 2. — Mass (solid line) and volume (dashed 
line) weighted neutral hydrogen fraction as a func- 
tion of redshift for the simulation A4 from Table 
1. Notice that the peak in the time derivative of 
the mean free path (vertical dashed line) closely 
corresponds to the moment when the mean mass 
(volume) weighted neutral fraction reaches a value 
of 10~ 2 (1CT 3 ) respectively. 

determined by (relatively) slowly evolving Lyman- 
limit systems (Miralda-Escude 2003a). 

It is, however, tempting to try to assign a value 
for the "redshift of reionization" . Such a value, 
in order not to be completely arbitrary, must be 
related to the physics of the reionization process. 
For example, it seems to be natural to identify the 
moment of reionization with the overlap stage, but 
even it takes place over a sizable redshift interval 
Az ~ 1 (Gnedin 2000), so we would need a bet- 
ter definition if we are to assign a value to "the 
redshift of reionization" . 

Fortunately, such a definition exists. Figure 1 
shows the mean free path to ionizing radiation and 
its time derivative as a function of redshift for the 
fiducial simulation (A4) with euv = 1 x 10~ 6 . The 
time derivative of the mean free path has a well- 
defined peak in the middle of the overlap stage, 
which is a natural moment to identify with "the 
redshift of reionization" zrei- 

However, this definition is not entirely practi- 
cal, since the time derivative of the mean free path 
cannot be observed directly. Instead, a more eas- 
ily observable quantity is the mean neutral frac- 
tion. Figure 2 shows the mean mass and volume 
weighted neutral hydrogen fraction for the fiducial 



simulation. The moment of reionization closely 
corresponds to the time when the mean mass (vol- 
ume) weighted neutral fraction reaches a value of 
10~ 2 (10~ 3 ) respectively - one can consider that 
moment as an alternative definition to "the red- 
shift of reionization" . 

The latter definition is more practical, but it 
is a subject to an important clause: while the 
volume weighted neutral fraction is reliably com- 
puted in the simulation, the mass weighted one 
depends on the numerical resolution. In particu- 
lar, the simulations presented in this paper do not 
resolve the damped Lyman-alpha systems, so the 
quoted above value of 1% for the mass weighted 
neutral fraction does not include the neutral gas 
locked in the damped Lyman-alpha systems. In 
this respect, the volume weighted number is more 
robust and should be used as a main definition of 
"the redshift of reionization" . 

4. Results 

4.1. Fitting the Mean Transmitted Flux 
Data 

Figure 3 shows the mean transmitted flux as a 
function of redshift for the three simulations with 
different values of the euv parameter from the set 
B4. The observational data shows with open cir- 
cles were obtained from White et al. (2003) by 
averaging the mean transmitted flux at a given 
redshift interval. The vertical error-bars are er- 
rors of the mean (not mean errors!). The horizon- 
tal error-bars are simply the width of the redshift 
interval over which the mean transmitted flux is 
computed. The last data point (without the ver- 
tical error-bar) is likely to be contaminated by a 
foreground galaxy and is not included in this anal- 
ysis. Filled squares show the data from Songaila 
(2004). In the latter case I do not show the error- 
bars for clarity, but they are comparable to that 
of White et al. (2003). 

Two gray lines are the solid black line simply 
shifted horizontally. As one can see, the change in 
the euv parameter simply translates into the shift 
in redshift. This is not surprising given that in the 
ACDM model clustering proceeds hierarchically, 
in a quasi self similar manner. In other words, 
one would expect that the dependence of the mean 
transmitted flux on the emissivity parameter euv 
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Fig. 3. — Mean transmitted flux as a function of 
redshift for simulations from the set B4. Solid, 
long-dashed, and short-dashed black lines corre- 
spond to simulations with euv = 2.5 x 10~ 6 , 
4.0 x 10~ 6 , and 2.0 x 10~ 6 respectively. Two gray 
lines are the solid black line simply shifted hor- 
izontally to match other two lines. Open circles 
with error-bars show the observational data from 
White et al. (2003) and filled squares are the data 
of Songaila (2004; the error-bars are not shown for 
these points for clarity). 

and redshift z enters, to the first order, only as 

<FH*) = (l + ^[e uv (l + (1) 

where g is a function of one argument, and a and 
(3 depend on the slope of the primordial power 
spectrum at the scale of interest, and may also 
depend on the details of the temperature-density 
relation. But because I am only concerned with a 
narrow redshift range around z — 6 (for example, 
z = 5.8-^6.5 in Fig. 3), a and (3 from equation (1) 
can be considered constant in the first approxima- 
tion. Fig. 3 indicates that for the range of redshifts 
considered (3 is small, less than about 0.5 in the ab- 
solute value. The specific reason for such a small 
value is not entirely clear, and may be a simple 
coincidence: if euv is reduced, the amount of ion- 
izing radiation in the post-overlap stage is also re- 
duced, but at the same time the post-overlap stage 
is delayed to lower redshifts, when densities are 
lower. The two effects appear to approximately 
cancel each other, leading to a small value of (3. 

The main result of Fig. 3 is that a change in the 
emissivity parameter euv by a factor of / simply 
translates into a rescale of (one plus) redshift by 
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Fig. 4. — Mean transmitted flux as a function of 
redshift for best-fit simulations from sets A4 (solid 
line), B4 (long-dashed line), and C4 (short-dashed 
line) . 

a factor f 1 ^. Note that such rescaling does not 
affect the vertical amplitude of the mean transmit- 
ted flux curve. As is well known, the mean absorp- 
tion after the post-overlap stage is dominated by 
the Lyman-limit systems (Miralda-Escudc 2003b) . 
Thus, in order to reproduce the observed run of 
the mean transmitted flux with redshift, a simu- 
lation must have (i) the right value for the emis- 
sivity parameter euv (i-e. the right value for the 
redshift of reionization) , to reproduce the drop- 
off at z ~ 6, and (ii) the correct abundance of 
the Lyman-limit systems. The latter is crucially 
dependent on the abundance of low mass virial- 
ized objects, and is not completely controlled by 
the parameter euv, but also depends on the mass 
function of virialized objects, i.e. on cosmological 
parameters. Therefore, if a simulation has wrong 
values of cosmological parameters, it might not 
be able to reproduce the amplitude of the mean 
transmitted flux as a function of redshift for any 
value of the single free parameter euv- 

This is illustrated in Figure 4 which shows the 
results from three simulation sets (A4, B4, and 
C4), each rescaled for the best value of the emis- 
sivity parameter euv (the simulation set C4 has 
not been continued beyond z — 5.5). As one can 
see, not all of the simulation sets are able to re- 
produce the amplitude of the mean transmitted 
flux curve. For example, the set C4 predicts a 
factor of 10 more absorption: in that cosmologi- 
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cal model (fi m = 0.35, n = 0.97) the amount of 
small-scale power is significantly larger than, say, 
in the WMAP model, and the mean distance be- 
tween Lyman-limit systems is a factor of 10 less 
than the observed value. Sets A4 and B4, on the 
other hand, provide a reasonable, although not 
perfect, fit to the data. In fact, it appears that the 
set A4 (the WMAP model) fits the data best for 
z > 5.2, but goes somewhat below the data points 
at lower redshifts. This is expected, since due to 
the limited size of the computational box, simula- 
tions should unpredict the abundance of ionizing 
sources at lower redshifts. The simulation set B4 
has a somewhat higher mean transmitted flux, al- 
though it is marginally consistent with the data 
given the cosmic variance. 

There is no easy way to reconcile model C4 
(and, perhaps, model B4) with the data. The 
abundance of the Lyman-limit systems depends 
not only on eyv, but also on the cosmological 
model, and, as long as euv an d cosmological pa- 
rameters arc specified, there is no free parameter 
left to adjust this abundance. The only way to 
make, say, model B4 to fit the data better would be 
to include additional Lyman-limit absorption pos- 
teriori, by hand - for example, by postulating an 
entirely different class of objects that contribute 
to the mean free path but do not form within the 
hierarchical clustering framework. For model C4, 
however, one would have to postulate that some 
of the Lyman-limit system that form in the simu- 
lation do not form in reality. 

4.2. Measuring the Redshift of Reioniza- 
tion 

Using the fact that the WMAP model fits the 
observed data, one can determine how much the 
redshift of reionization can be shifted while still 
preserving the acceptable fit to the data. Figure 5 
shows the range of values for the redshift of reion- 
ization for models A4 and B4 (although the lat- 
ter is admittedly a worse fit to the data) obtained 
by adding linearly two different components: (i) 
an average observational error-bar in z-dircction 
(0.085), and (ii) the average uncertainty to the 
best fit value of zrei for the fiducial model. With 
these three uncertainties included, a measurement 
of the redshift of reionization can be made: 

zrei = 6.1 ±0.3 (95% CL). (2) 
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Fig. 5. — Mean transmitted flux as a function 
of redshift for the best-fit simulations from sets 
A4 (dark gray) and B4 (medium gray). Shaded 
regions illustrate the uncertainty due to observa- 
tional errors. Arrows with error-bars give the red- 
shift of reionization plus its la uncertainty for each 
of the simulation sets. 

If other cosmological models are also considered 
(for example, model B4), then an additional "sys- 
tematic" uncertainty is introduced due to small 
differences in the redshifts of reionization of differ- 
ent cosmological models. This uncertainty, how- 
ever, is rather small, as can be seen from Fig. 5, 
and is ignored here. 

4.3. Sanity Checks 

Before I can proceed further, several sanity 
tests must be completed. The first two are shown 
in Figure 6. The light gray line is again the best- 
fit fiducial WMAP model, with error-bars showing 
the rms fluctuation in the mean transmitted flux 
(error of the mean) - this is consistent with (al- 
though slightly smaller than) the observed values. 
This is not surprising, as a finite (and relatively 
small) size of the computational volume will lead 
to an underestimate of the cosmic variance. 

How significant this overestimate? In order to 
reduce the cosmic variance at the box size scale, 
the specific realizations of initial conditions are se- 
lected to reproduce the correct power spectrum 
of initial perturbations on the fundamental mode 
(Gnedin & Hamilton 2002), which mitigates the 
uncertainty due to cosmic variance by about a fac- 
tor of 10 or more. To demonstrate that numerical 
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Fig. 6. — Tests of numerical convergence and cos- 
mic variance. Light gray line with error-bars is the 
best fit A4 set with the expected rms fluctuation in 
the mean transmitted flux over a redshift interval 
Az = 0.17 (which is the the average redshift in- 
terval over which the transmitted flux is averaged 
in the observational data). The dark gray long- 
dashed line is the best-fit mean transmitted flux 
for the A8 set (twice larger simulation volume). 

convergence is achieved for the mean transmitted 
flux, and that cosmic variance does not affect my 
results, I show with the dark gray line in Fig. 6 the 
best-fit simulation from the A8 set (the same cos- 
mological model, twice larger box). The fact that 
the two simulations with the two box sizes are suf- 
ficiently close argues for the numerical convergence 
of the simulation result, although, admittedly, this 
is not a 100% rigorous test (three different reso- 
lutions would be required for a complete test, but 
this is currently beyond the existing computer ca- 
pabilities). 

The recent results from the WMAP measure- 
ment of the cross-correlation between the temper- 
ature and polarization anisotropies of the CMB 
indicate the large value of the Compton optical 
depth t = 0.12 ±0.06 (Kogut et al. 2003, Tegmark 
et al. 2004). Does this observation rule out reion- 
ization at z ~ 6? Figure 7 shows the ioniza- 
tion histories for the fiducial model (which has 
r = 0.06) together with three other posteriori ion- 
ization histories which all have r = 0.2. Because 
these histories are adjusted to match the fiducial 
model before the volume averaged ionized fraction 
approaches 99.9%, they all have identical redshifts 
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Fig. 7. — Ionization histories (volume averaged 
ionized fraction as a function of redshift) for the 
fiducial model (black line), and three arbitrary 
ionization histories (grey lines) with r = 0.2. All 
four model have the identical ionization redshift of 
zrei =6.1. 

of reionization and identical mean transmitted flux 
curves at z < 6.1. Thus, the WMAP measure- 
ment by itself neither constraints the redshift of 
reionization nor contradicts the SLOAN data, but 
rather indicates that the early ionization history 
might have been quite complex. 

5. Conclusions 

I have shown that advanced self-consistent cos- 
mological simulation with radiative transfer are 
able to reproduce the observed data on the evolu- 
tion of the mean transmitted flux from the SLOAN 
survey when the values of cosmological parame- 
ters that best fit the WMAP data are adopted for 
the underlying cosmology. Thus, a consistent pic- 
ture is emerging in which the universe reionizes 
at 6.1 ± 0.3 after possibly a prolonged period of 
partial reionization. 

I am grateful to the anonymous referee for cor- 
recting errors in the original manuscript, and for 
other valuable comments. This work was sup- 
ported by NSF grant AST-0134373 and by Na- 
tional Computational Science Alliance under grant 
AST-020018N and utilized SGI Origin 2000 array 
and IBM P690 array at the National Center for 
Supcrcomputing Applications. 
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